In [1]:
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import minimize 
from datetime import datetime as dt
from eod import EodHistoricalData
from nelson_siegel_svensson import NelsonSiegelSvenssonCurve
from nelson_siegel_svensson.calibrate import calibrate_nss_ols, calibrate_ns_ols

import matplotlib

# import matplotlib.pyplot as plt
# from mpl_toolkits.mplot3d import axes3d  

from matplotlib.backends.backend_pgf import FigureCanvasPgf
matplotlib.backend_bases.register_backend('pdf', FigureCanvasPgf)

matplotlib.rcParams.update({
    "pgf.texsystem": "pdflatex",
    'font.family': 'serif',
    'text.usetex': True,
    'font.size': 14,
    'legend.fontsize': 13,
    'pgf.rcfonts': False
})
In [2]:
from scipy.stats import norm
import scipy.integrate as integrate
from scipy.special import iv

def blackScholes(r, S, K, T, sigma, type="c"):
    d1 = (np.log(S/K) + (r + sigma**2/2)*T)/(sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    if type == "c":
        price = S*norm.cdf(d1, 0, 1) - K*np.exp(-r*T)*norm.cdf(d2, 0, 1)
    elif type == "p":
        price = K*np.exp(-r*T)*norm.cdf(-d2, 0, 1) - S*norm.cdf(-d1, 0, 1)
    return price

def switchingCall(X0,K,T,r,s0,s1,l0,l1,pi=-1):
    
    if pi == -1:
        pi = l1/(l0+l1)
    
    def s(x, T):
        return s0**2 * x + s1**2 * (T-x)

    def g0(x, T):
        return np.sqrt(l0*l1*x/(T-x))

    def g1(x, T):
        return np.sqrt(l0*l1*(T-x)/x)

    def h(x, T):
        return np.sqrt(l0*l1*x*(T-x))

    def f0(x, T):
        return np.exp(-l0 * x - l1 * (T-x))*(g0(x,T)*iv(1, 2*h(x,T)) + l0 * iv(0, 2*h(x,T)))

    def f1(x,T):
        return np.exp(-l0 * x - l1 * (T-x))*(g1(x,T)*iv(1, 2*h(x,T)) + l1 * iv(0, 2*h(x,T)))

    def integrand0(x, K, T):
        return blackScholes(r, X0, K, T, np.sqrt(s(x,T)/T)) * f0(x, T)

    def integrand1(x, K, T):
        return blackScholes(r, X0, K, T, np.sqrt(s(x,T)/T)) * f1(x, T)

    def C0(K, T):
        return blackScholes(r, X0, K, T, s0)*np.exp(-l0*T) + integrate.quad(integrand0, 0, T, args=(K,T))[0]

    def C1(K, T):
        return blackScholes(r, X0, K, T, s1)*np.exp(-l0*T) + integrate.quad(integrand1, 0, T, args=(K,T))[0]
    
#     print(C0(K,T), C1(K,T))
    return pi * C0(K,T) + (1-pi) * C1(K,T)
In [3]:
yield_maturities = np.array([1/12, 2/12, 3/12, 4/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30])
yields = np.array([5.51,5.56,5.56,5.60,5.49,5.43,5.02,4.72,4.45,4.41,4.33,4.59,4.42]).astype(float)/100
In [4]:
#NSS model calibrate
curve_fit, status = calibrate_nss_ols(yield_maturities,yields) 
curve_fit
Out[4]:
NelsonSiegelSvenssonCurve(beta0=0.048410223860134056, beta1=0.008319582817567225, beta2=-0.001366128289200936, beta3=-0.021638282142028585, tau1=2.0, tau2=5.0)
In [5]:
t = np.linspace(0, 30, 100)
plt.plot(t, curve_fit(t)*100, color='black', label='Nelson-Siegel-Svensson Model')
plt.scatter(yield_maturities, yields*100, facecolors='none', edgecolors='black', label='Treasury Par Yield Curve Rate')
plt.grid(linewidth = 0.5)
plt.legend()
plt.ylabel("Rate (\%)")
plt.xlabel("Tenor (Years)")

# plt.savefig("yield.pgf")
In [6]:
curve_fit, status = calibrate_ns_ols(yield_maturities,yields) 
print(curve_fit)

t = np.linspace(0, 30, 100)
plt.plot(t, curve_fit(t)*100, color='black', label='Nelson-Siegel Model')
plt.scatter(yield_maturities, yields*100, facecolors='none', edgecolors='black', label='Treasury Par Yield Curve Rate')
plt.grid(linewidth = 0.5)
plt.legend()
plt.ylabel("Rate (\%)")
plt.xlabel("Tenor (Years)")

# plt.savefig("yield2.pgf")
NelsonSiegelCurve(beta0=0.04410353509255659, beta1=0.013119399209004352, beta2=-0.011705340534500869, tau=2.0)
In [7]:
import yfinance as yf

def option_chains(ticker):
    
    asset = yf.Ticker(ticker)
#     print(asset)
    expirations = asset.options
    
#     print(expirations[:10])
    
    market_prices = {}
    
    for expiration in expirations:
        opt = asset.option_chain(expiration)
        chain = opt.calls
    
#         print(chain)
        
        expiration = pd.to_datetime(expiration) + pd.DateOffset(hours=23, minutes=59, seconds=59)
        
        market_prices[expiration] = {}
        market_prices[expiration]['strike'] = chain['strike'].tolist()
#         market_prices[expiration]['price'] = (chain['bid']/2 + chain['ask']/2).tolist()
        market_prices[expiration]['price'] = chain['lastPrice'].tolist()

    
    return market_prices

options = option_chains("SPY")

# fig = plt.figure()
# ax = fig.add_subplot(projection='3d')
# for i,v in options.items():
#     ax.scatter(v['maturity'], v['strike'], v['lastPrice'])

all_strikes = [v['strike'] for i,v in options.items()]
common_strikes = set.intersection(*map(set,all_strikes))
# common_strikes.update([420,430,480,490])
print('Number of common strikes:', len(common_strikes))
common_strikes = sorted(common_strikes)

prices = []
maturities = []
for date, v in options.items():
    maturities.append((date - dt.today()).days/365.25)
    price = [v['price'][i] for i,x in enumerate(v['strike']) if x in common_strikes]
    prices.append(price)
price_arr = np.array(prices, dtype=object)
np.shape(price_arr)
volSurface = pd.DataFrame(price_arr, index = maturities, columns = common_strikes)
volSurface = volSurface.iloc[(volSurface.index > 0.04) & (volSurface.index <= 1)]
volSurface
Number of common strikes: 6
Out[7]:
415.0 420.0 425.0 430.0 440.0 450.0
0.052019 16.37 12.36 8.84 5.8 1.74 0.27
0.071184 17.57 13.23 9.83 6.99 2.57 0.6
0.090349 18.55 14.45 11.29 8.24 3.55 1.06
0.109514 19.21 15.43 12.19 8.75 4.21 1.57
0.128679 20.6 16.76 13.27 10.12 5.09 2.02
0.205339 23.83 20.03 16.66 13.22 7.9 4
0.243669 24.38 20.64 17.31 14.29 8.5 4.65
0.301164 26.5 23.16 19.62 16.4 10.66 6.25
0.454483 31.14 28.47 25.62 22.13 16.05 10.85
0.490075 33.46 31 26 21.6 16.36 11.9
0.722793 42.24 40.05 34.85 30.55 24.21 19.02
0.741958 42.95 39.87 35.7 33.26 25.23 19
0.971937 47.75 46.74 40.5 39.33 31.77 27.06
In [8]:
def convert_prices_to_chain(prices):
    output = pd.DataFrame()
    for expiry, pair in prices.items():
#         print(expiry, pair)
        for strike, price in zip(pair['strike'], pair['price']):
            output.loc[((expiry - dt.today()).days/365.25), strike] = price
            
    output = output.reindex(sorted(output.columns), axis=1)
    
    return output

volSurface = convert_prices_to_chain(options)
volSurface = volSurface.loc[:, (volSurface.columns >= 400.0) & (volSurface.columns <= 500.0)]
volSurface = volSurface.iloc[(volSurface.index > 0.2)]
volSurface
Out[8]:
400.0 401.0 402.0 403.0 404.0 405.0 406.0 407.0 408.0 409.0 ... 491.0 492.0 493.0 494.0 495.0 496.0 497.0 498.0 499.0 500.0
0.205339 35.96 34.72 34.63 31.78 32.49 32.32 31.11 30.29 30.39 29.28 ... 0.14 0.13 0.15 0.11 0.11 0.12 0.11 0.11 0.11 0.09
0.243669 37.92 NaN NaN NaN NaN 32.50 NaN NaN 30.21 28.84 ... 0.20 0.18 0.15 0.18 0.15 NaN 0.14 0.13 0.12 0.11
0.301164 37.85 NaN NaN NaN NaN 34.71 NaN NaN NaN NaN ... NaN NaN NaN NaN 0.27 NaN NaN NaN NaN 0.19
0.454483 43.66 NaN NaN NaN NaN 40.35 NaN NaN NaN NaN ... NaN NaN NaN NaN 0.92 NaN NaN NaN NaN 0.70
0.490075 44.68 NaN NaN NaN NaN 40.29 NaN NaN NaN NaN ... NaN NaN NaN NaN 1.10 NaN NaN NaN NaN 0.85
0.722793 53.26 NaN NaN NaN NaN 47.08 NaN NaN NaN NaN ... NaN NaN NaN NaN 3.59 NaN NaN NaN NaN 2.90
0.741958 53.60 51.96 50.57 50.48 58.63 57.76 56.72 53.84 55.16 55.00 ... NaN NaN NaN NaN 4.05 NaN NaN NaN NaN 2.85
0.971937 61.94 NaN NaN NaN NaN 54.68 NaN NaN NaN NaN ... NaN NaN NaN NaN 7.75 NaN NaN NaN NaN 6.14
1.221081 64.68 NaN NaN NaN NaN 61.74 NaN NaN NaN NaN ... NaN NaN NaN NaN 11.47 NaN NaN NaN NaN 10.58
1.297741 64.93 NaN NaN NaN NaN 62.70 NaN NaN NaN NaN ... NaN NaN NaN NaN 13.22 NaN NaN NaN NaN 10.06
1.470226 85.10 NaN NaN NaN NaN 67.44 NaN NaN NaN NaN ... NaN NaN NaN NaN 16.41 NaN NaN NaN NaN 16.00
1.719370 75.00 NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN 22.12 NaN NaN NaN NaN 21.57
2.217659 83.90 NaN NaN NaN NaN 77.00 NaN NaN NaN NaN ... NaN NaN NaN NaN 29.36 NaN NaN NaN NaN 30.00
2.294319 83.87 NaN NaN NaN NaN 78.00 NaN NaN NaN NaN ... NaN NaN NaN NaN 32.64 NaN NaN NaN NaN 29.90

14 rows × 110 columns

In [9]:
print(volSurface.columns.values)
[400.  401.  402.  403.  404.  405.  406.  407.  408.  409.  410.  411.
 412.  413.  414.  415.  416.  417.  418.  419.  420.  421.  422.  422.5
 423.  424.  425.  426.  427.  427.5 428.  429.  430.  431.  432.  432.5
 433.  434.  435.  436.  437.  437.5 438.  439.  440.  441.  442.  442.5
 443.  444.  445.  446.  447.  447.5 448.  449.  450.  451.  452.  452.5
 453.  454.  455.  456.  457.  457.5 458.  459.  460.  461.  462.  462.5
 463.  464.  465.  466.  467.  468.  469.  470.  471.  472.  473.  474.
 475.  476.  477.  478.  479.  480.  481.  482.  483.  484.  485.  486.
 487.  488.  489.  490.  491.  492.  493.  494.  495.  496.  497.  498.
 499.  500. ]
In [10]:
S0 = 445

from scipy.stats import norm
N = norm.cdf

def bs_call(S, K, T, r, vol):
    d1 = (np.log(S/K) + (r + 0.5*vol**2)*T) / (vol*np.sqrt(T))
    d2 = d1 - vol * np.sqrt(T)
    return S * norm.cdf(d1) - np.exp(-r * T) * K * norm.cdf(d2)

def bs_vega(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    return S * norm.pdf(d1) * np.sqrt(T)

def find_vol(target_value, S, K, T, r, *args):
    MAX_ITERATIONS = 200
    PRECISION = 1.0e-5
    sigma = 0.5
    for i in range(0, MAX_ITERATIONS):
        price = bs_call(S, K, T, r, sigma)
        vega = bs_vega(S, K, T, r, sigma)
        diff = target_value - price  # our root
        if (abs(diff) < PRECISION).all() :
            return sigma
        sigma = sigma + diff/vega # f(x) / f'(x)
    return sigma # value wasn't found, return best guess so far

vect = np.vectorize(find_vol)
In [11]:
# Convert our vol surface to dataframe for each option price with parameters
volSurfaceLong = volSurface.melt(ignore_index=False).dropna().reset_index()
volSurfaceLong.columns = ['maturity', 'strike', 'price']
# Calculate the risk free rate for each maturity using the fitted yield curve
volSurfaceLong['rate'] = volSurfaceLong['maturity'].apply(curve_fit)
# volSurfaceLong['marketIV'] = volSurfaceLong.apply(lambda x: vect(x.price, S0, x.strike, x.maturity, x.rate), axis=1)
volSurfaceLong
Out[11]:
maturity strike price rate
0 0.205339 400.0 35.96 0.056011
1 0.243669 400.0 37.92 0.055798
2 0.301164 400.0 37.85 0.055485
3 0.454483 400.0 43.66 0.054694
4 0.490075 400.0 44.68 0.054519
... ... ... ... ...
566 1.297741 500.0 10.06 0.051261
567 1.470226 500.0 16.00 0.050717
568 1.719370 500.0 21.57 0.050007
569 2.217659 500.0 30.00 0.048820
570 2.294319 500.0 29.90 0.048662

571 rows × 4 columns

In [12]:
S0 = 428.8
r = volSurfaceLong['rate'].to_numpy('float')
K = volSurfaceLong['strike'].to_numpy('float')
tau = volSurfaceLong['maturity'].to_numpy('float')
P = volSurfaceLong['price'].to_numpy('float')
In [13]:
params = {"sig0" : {"x0": 0.16, "lbub": [1e-2,1]},
          "sig1" : {"x0": 0.07, "lbub": [1e-2,1]},
          "l0"   : {"x0": 0.76, "lbub": [1e-2,5]},
          "l1"   : {"x0": 1.97, "lbub": [1e-2,5]}
          }

x0 = [param["x0"] for key, param in params.items()]
bnds = [param["lbub"] for key, param in params.items()]

# switchingCall(X0,K,T,r,s0,s1,l0,l1)

def callback(xk):
    print(f"Iteration: {callback.nfev}")
    print(f"Current solution: {xk}")
    callback.nfev += 1

callback.nfev = 0

def SqErr(x):
    

        
    sig0, sig1, l0, l1 = [param for param in x]
    
    err = np.sum([ (P_i-switchingCall(S0, K_i, tau_i, r_i, sig0, sig1, l0, l1))**2 /len(P) \
                  for P_i, K_i, tau_i, r_i in zip(P, K, tau, r)])
    
    pen = np.sum( [(x_i-x0_i)**2 for x_i, x0_i in zip(x, x0)] )
          
    return err + pen

print("About to calibrate")

result = minimize(SqErr, x0, tol = 1e-5, method='L-BFGS-B', options={'maxiter': 1e4 }, bounds=bnds, callback=callback)
sig0, sig1, l0, l1 = [param for param in result.x]
sig0, sig1, l0, l1

# (0.10074434838728742, 0.15631085210777004, 0.807983501745533, 2.0304931134722355)

# About to calibrate
# Iteration: 0
# Current solution: [0.01 0.01 5.   0.01]
# Iteration: 1
# Current solution: [0.01       0.10952676 2.60712973 1.35982573]
# Iteration: 2
# Current solution: [0.01       0.19772305 2.45767573 1.42899444]
# Iteration: 3
# Current solution: [0.01       0.2280412  2.34995689 1.48091238]
# Iteration: 4
# Current solution: [0.01       0.22237417 2.29474563 1.51001524]
# Iteration: 5
# Current solution: [0.01       0.20246025 1.83737575 1.74714761]
# Iteration: 6
# Current solution: [0.01111231 0.20483723 1.63835345 1.84657843]
# Iteration: 7
# Current solution: [0.03534615 0.21478326 1.38540795 1.951807  ]
# Iteration: 8
# Current solution: [0.04896322 0.21736113 1.17978456 2.043398  ]
# Iteration: 9
# Current solution: [0.05589635 0.2177989  0.94020518 2.1540299 ]
# Iteration: 10
# Current solution: [0.06869119 0.22015591 0.79607144 2.21474713]
# Iteration: 11
# Current solution: [0.07962909 0.22192583 0.63965716 2.28260622]
# Iteration: 12
# Current solution: [0.08410475 0.22468342 0.55137829 2.29334664]
# Iteration: 13
# Current solution: [0.08545409 0.22915719 0.5144569  2.24198814]
# Iteration: 14
# Current solution: [0.08127998 0.24277412 0.43500788 2.06614759]
# Iteration: 15
# Current solution: [0.07894165 0.24418271 0.45237824 2.02970636]
# Iteration: 16
# Current solution: [0.0775697  0.24716911 0.43740063 1.98573046]
# sig0, sig1, l0, l1 = (0.07756969778166047, 0.247169113443559, 0.4374006282788692, 1.9857304635023587)
About to calibrate
Iteration: 0
Current solution: [0.16383145 0.07424196 0.75657906 1.97820498]
Iteration: 1
Current solution: [0.16592249 0.07659084 0.75460704 1.98277015]
Iteration: 2
Current solution: [0.16587679 0.07652628 0.75456612 1.98270935]
Iteration: 3
Current solution: [0.16580624 0.07639306 0.75428318 1.98271817]
Iteration: 4
Current solution: [0.1656938  0.07609079 0.75324131 1.98300832]
Iteration: 5
Current solution: [0.16552656 0.07540177 0.75011902 1.98417411]
Iteration: 6
Current solution: [0.1652965  0.07379826 0.74153393 1.98777724]
Iteration: 7
Current solution: [0.16504642 0.07017283 0.72000551 1.99735884]
Iteration: 8
Current solution: [0.16496461 0.06272966 0.67321163 2.0188374 ]
Iteration: 9
Current solution: [0.16605768 0.05085295 0.61232799 2.04509635]
Iteration: 10
Current solution: [0.16616216 0.0530357  0.62433756 2.03980222]
Iteration: 11
Current solution: [0.16632257 0.05630798 0.6367544  2.03412554]
Iteration: 12
Current solution: [0.1662957  0.0572916  0.6373558  2.03383143]
Iteration: 13
Current solution: [0.1660157  0.05901155 0.6348142  2.03505249]
Iteration: 14
Current solution: [0.16512717 0.06265337 0.62729144 2.03877971]
Iteration: 15
Current solution: [0.16364385 0.06728713 0.61473343 2.04501553]
Iteration: 16
Current solution: [0.16109842 0.07935921 0.62681562 2.03703891]
Iteration: 17
Current solution: [0.15894926 0.08286976 0.60337518 2.04324121]
Iteration: 18
Current solution: [0.15900836 0.08251891 0.60475824 2.04139208]
Iteration: 19
Current solution: [0.15901135 0.08234573 0.60720606 2.03612614]
Iteration: 20
Current solution: [0.15874187 0.08310705 0.60874033 2.02597484]
Iteration: 21
Current solution: [0.15747949 0.08728243 0.60738006 1.99869935]
Iteration: 22
Current solution: [0.15150136 0.10761231 0.59070681 1.91466738]
Iteration: 23
Current solution: [0.14878163 0.11621042 0.57911009 1.88478659]
Iteration: 24
Current solution: [0.14281347 0.13132854 0.5283959  1.81875637]
Iteration: 25
Current solution: [0.14350182 0.12916196 0.53063472 1.82675944]
Iteration: 26
Current solution: [0.143944   0.12815811 0.53598363 1.83112665]
Iteration: 27
Current solution: [0.14314842 0.13017471 0.53558182 1.82084361]
Iteration: 28
Current solution: [0.14047352 0.13712548 0.53138074 1.80044245]
Iteration: 29
Current solution: [0.13732007 0.14526308 0.52628622 1.77853312]
Iteration: 30
Current solution: [0.13504951 0.15086456 0.52319791 1.76463436]
Iteration: 31
Current solution: [0.12821123 0.16693966 0.51528946 1.73065036]
Iteration: 32
Current solution: [0.12650539 0.16983148 0.51559311 1.73375405]
Iteration: 33
Current solution: [0.12041979 0.18255591 0.51182894 1.72622768]
Iteration: 34
Current solution: [0.1024681  0.21987735 0.49550008 1.78391922]
Iteration: 35
Current solution: [0.11272122 0.19701024 0.52361094 1.83135075]
Iteration: 36
Current solution: [0.11249968 0.19831077 0.53425235 1.79785939]
Iteration: 37
Current solution: [0.10900122 0.20553972 0.52173099 1.80857967]
Iteration: 38
Current solution: [0.10866632 0.20618749 0.51931773 1.80659175]
Iteration: 39
Current solution: [0.1087006  0.2060978  0.51864845 1.80368293]
Out[13]:
(0.10870059646871238,
 0.20609779580245238,
 0.5186484518816581,
 1.8036829298762445)
In [38]:
# sig0, sig1, l0, l1 = (0.07756969778166047, 0.247169113443559, 0.4374006282788692, 1.9857304635023587)
# sig0, sig1, l0, l1 = (0.07580256, 0.25515649, 0.54693518, 1.91751383)
In [39]:
vectorized_call = np.vectorize(switchingCall)

switch_prices = vectorized_call(S0, K, tau, r, sig0, sig1, l0, l1)
volSurfaceLong['switch_price'] = switch_prices
In [40]:
# import pickle

# with open('best_calibration.pkl', 'wb') as f:
#     pickle.dump(volSurfaceLong, f)
    
# with open('best_calibration_final.pkl', 'rb') as f:
#     volSurfaceLong = pickle.load(f)
In [41]:
import plotly.graph_objects as go
from plotly.graph_objs import Surface
from plotly.offline import iplot, init_notebook_mode
init_notebook_mode()

x = np.linspace(400, 500, 100)
y = np.linspace(0.1, 2.3, 23)

X, Y = np.meshgrid(x, y)

Z = vectorized_call(S0, X, Y, curve_fit(Y), sig0, sig1, l0, l1)
In [42]:
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.add_scatter3d(x=volSurfaceLong.strike, y=volSurfaceLong.maturity, z=volSurfaceLong.price, mode='markers')
fig.update_layout(
    title_text='Market Prices (Mesh) vs Calibrated Regime-switching Prices (Markers)',
    scene = dict(yaxis_title='TIME (Years)',
                    xaxis_title='STRIKES (Pts)',
                    zaxis_title='INDEX OPTION PRICE (Pts)'),
    height=800,
    width=800
)
fig.show()
In [43]:
# print(list(volSurface.index.values), list(volSurface.columns.values))
In [44]:
x = np.linspace(400, 500, 50)
y = np.linspace(0.005, 2.3, 25)

X, Y = np.meshgrid(x, y)
# print(X.size,Y.size)
Z = vectorized_call(S0, X, Y, curve_fit(Y), sig0, sig1, l0, l1)
In [45]:
X.shape, Y.shape, Z.shape
Out[45]:
((25, 50), (25, 50), (25, 50))
In [46]:
Z[0]
Out[46]:
array([2.89545395e+01, 2.69114967e+01, 2.48685045e+01, 2.28256509e+01,
       2.07831460e+01, 1.87414458e+01, 1.67014529e+01, 1.46647926e+01,
       1.26341267e+01, 1.06134090e+01, 8.60799576e+00, 6.62567525e+00,
       4.69010226e+00, 2.90164587e+00, 1.49780953e+00, 6.73690268e-01,
       3.20146557e-01, 1.80057990e-01, 1.07679153e-01, 6.26563634e-02,
       3.47288691e-02, 1.82769573e-02, 9.12336579e-03, 4.31660452e-03,
       1.93487395e-03, 8.21388381e-04, 3.30183665e-04, 1.25676182e-04,
       4.52965269e-05, 1.54615559e-05, 4.99932261e-06, 1.53164586e-06,
       4.44770719e-07, 1.22462794e-07, 3.19842961e-08, 7.92722071e-09,
       1.86532109e-09, 4.16907748e-10, 8.85506796e-11, 1.78824304e-11,
       3.43530403e-12, 6.28105003e-13, 1.09359003e-13, 1.81409673e-14,
       2.86866194e-15, 4.32654524e-16, 6.22695961e-17, 8.55684030e-18,
       1.12326544e-18, 1.40932465e-19])
In [47]:
plt.figure(figsize=(8, 6), dpi=80)
ax = plt.axes(projection='3d')
ax.plot_wireframe(X, Y, Z, color='black', cstride=5, rstride=2, label='Calibrated Price Surface')
ax.scatter(volSurfaceLong.strike, volSurfaceLong.maturity, volSurfaceLong.price, label='Market Prices')
ax.set_xlabel('Strike (K)')
ax.set_ylabel('Time to Maturity (T)')
ax.set_zlabel('Price')
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.legend()
# plt.savefig("fittedpricesurface.pgf",bbox_inches='tight')
Out[47]:
<matplotlib.legend.Legend at 0x7f8818126dd8>
In [48]:
volSurfaceLong
Out[48]:
maturity strike price rate switch_price implVol marketIV error IVerror
0 0.205339 400.0 35.96 0.056011 36.317754 0.182439 0.174492 0.994867 0.045543
1 0.243669 400.0 37.92 0.055798 37.703242 0.181256 0.185341 -0.571620 -0.022040
2 0.301164 400.0 37.85 0.055485 39.736194 0.180502 0.147125 4.983341 0.226863
3 0.454483 400.0 43.66 0.054694 44.869682 0.180643 0.165353 2.770687 0.092470
4 0.490075 400.0 44.68 0.054519 46.000349 0.180788 0.164876 2.955123 0.096508
... ... ... ... ... ... ... ... ... ...
566 1.297741 500.0 10.06 0.051261 12.653258 0.137335 0.122129 25.777910 0.124504
567 1.470226 500.0 16.00 0.050717 15.076666 0.136410 0.141205 -5.770835 -0.033958
568 1.719370 500.0 21.57 0.050007 18.624106 0.135354 0.148976 -13.657366 -0.091433
569 2.217659 500.0 30.00 0.048820 25.749684 0.133735 0.150519 -14.167720 -0.111509
570 2.294319 500.0 29.90 0.048662 26.842108 0.133519 0.145372 -10.227065 -0.081536

571 rows × 9 columns

In [49]:
volSurfaceLong['implVol'] = vect(volSurfaceLong['switch_price'], S0, volSurfaceLong['strike'], \
                                    volSurfaceLong['maturity'], volSurfaceLong['rate'])

volSurfaceLong['marketIV'] = vect(volSurfaceLong['price'], S0, volSurfaceLong['strike'], \
                                    volSurfaceLong['maturity'], volSurfaceLong['rate'])

volSurfaceLong
Out[49]:
maturity strike price rate switch_price implVol marketIV error IVerror
0 0.205339 400.0 35.96 0.056011 36.317754 0.182439 0.174492 0.994867 0.045543
1 0.243669 400.0 37.92 0.055798 37.703242 0.181256 0.185341 -0.571620 -0.022040
2 0.301164 400.0 37.85 0.055485 39.736194 0.180502 0.147125 4.983341 0.226863
3 0.454483 400.0 43.66 0.054694 44.869682 0.180643 0.165353 2.770687 0.092470
4 0.490075 400.0 44.68 0.054519 46.000349 0.180788 0.164876 2.955123 0.096508
... ... ... ... ... ... ... ... ... ...
566 1.297741 500.0 10.06 0.051261 12.653258 0.137335 0.122129 25.777910 0.124504
567 1.470226 500.0 16.00 0.050717 15.076666 0.136410 0.141205 -5.770835 -0.033958
568 1.719370 500.0 21.57 0.050007 18.624106 0.135354 0.148976 -13.657366 -0.091433
569 2.217659 500.0 30.00 0.048820 25.749684 0.133735 0.150519 -14.167720 -0.111509
570 2.294319 500.0 29.90 0.048662 26.842108 0.133519 0.145372 -10.227065 -0.081536

571 rows × 9 columns

In [50]:
find_vol(45.38, 445,400,0.0566,0.01369)
Out[50]:
0.19557222628611118
In [51]:
x = np.linspace(400, 500, 50)
y = np.linspace(0.1, 2.3, 25)

X,Y = np.meshgrid(x,y)

Z = vect(vectorized_call(S0, X, Y, curve_fit(Y), sig0, sig1, l0, l1), S0, X, Y, curve_fit(Y))
In [52]:
fig2 = go.Figure(data=[go.Surface(x=X, y=Y, z=Z, opacity=0.55)])
fig2.update_layout(
    title_text='Implied Volatility',
    scene = dict(xaxis_title='TIME (Years)',
                    yaxis_title='STRIKES (Pts)',
                    zaxis_title='Implied Volatility'),
    height=800,
    width=800
)
fig2.show()
In [53]:
plt.figure(figsize=(8, 6), dpi=80)
ax = plt.axes(projection='3d')
ax.plot_wireframe(X, Y, Z, color='black', cstride=5, rstride=2)
# ax.scatter(volSurfaceLong.strike, volSurfaceLong.maturity, volSurfaceLong.price)
ax.set_xlabel('Strike (K)')
ax.set_ylabel('Time to Maturity (T)')
ax.set_zlabel('Implied Volatility')
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
# plt.savefig("fittedvolsurface.pgf", bbox_inches='tight')
In [54]:
volSurfaceLong['error'] = (volSurfaceLong['switch_price']-volSurfaceLong['price'])/volSurfaceLong['price']
In [55]:
volSurfaceLong
Out[55]:
maturity strike price rate switch_price implVol marketIV error IVerror
0 0.205339 400.0 35.96 0.056011 36.317754 0.182439 0.174492 0.009949 0.045543
1 0.243669 400.0 37.92 0.055798 37.703242 0.181256 0.185341 -0.005716 -0.022040
2 0.301164 400.0 37.85 0.055485 39.736194 0.180502 0.147125 0.049833 0.226863
3 0.454483 400.0 43.66 0.054694 44.869682 0.180643 0.165353 0.027707 0.092470
4 0.490075 400.0 44.68 0.054519 46.000349 0.180788 0.164876 0.029551 0.096508
... ... ... ... ... ... ... ... ... ...
566 1.297741 500.0 10.06 0.051261 12.653258 0.137335 0.122129 0.257779 0.124504
567 1.470226 500.0 16.00 0.050717 15.076666 0.136410 0.141205 -0.057708 -0.033958
568 1.719370 500.0 21.57 0.050007 18.624106 0.135354 0.148976 -0.136574 -0.091433
569 2.217659 500.0 30.00 0.048820 25.749684 0.133735 0.150519 -0.141677 -0.111509
570 2.294319 500.0 29.90 0.048662 26.842108 0.133519 0.145372 -0.102271 -0.081536

571 rows × 9 columns

In [56]:
volSurfaceLong[volSurfaceLong.error > 100]
Out[56]:
maturity strike price rate switch_price implVol marketIV error IVerror
In [57]:
volSurfaceLong['IVerror'] = (volSurfaceLong['implVol']-volSurfaceLong['marketIV'])/volSurfaceLong['marketIV']
In [58]:
volSurfaceLong[volSurfaceLong.IVerror > 20]
Out[58]:
maturity strike price rate switch_price implVol marketIV error IVerror
In [59]:
import plotly.express as px

fig = px.scatter_3d(volSurfaceLong, x='maturity', y='strike', z='error')
fig.show()
In [60]:
plt.figure(figsize=(8, 6), dpi=80)
ax = plt.axes(projection='3d')
ax.scatter(volSurfaceLong.strike, volSurfaceLong.maturity, volSurfaceLong.error)
ax.set_xlabel('Strike (K)')
ax.set_ylabel('Time to Maturity (T)')
ax.set_zlabel('Error (\%)')
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.view_init(30, 120)
# plt.savefig("calibrationerror.pgf", bbox_inches='tight')
In [63]:
import plotly
plotly.offline.init_notebook_mode()